## iterative model specification search ##
## load data ##

	rm(list = ls(all=TRUE))
	library(mnormt)
	library(lme4)
	library(stargazer)
	data <- read.table('fullData.txt', header = TRUE, sep = '|')
	data$district <- paste(data$state, data$district)

## main models ##
## simple measure (main text table 1) ##
	winsModel <- lmer((majvote/100) ~ lagDv + winsParty + wins + spendGapIMP + majControl + majInc + minControl + minInc + majorityQuality + minorityQuality + unopposed + majPresVote*presYear*majPresAgree + (1 | year) + (1 | district), data = data)
	summary(winsModel)

## four-part measure (appendix) ##
	fourPart <- lmer((majvote/100) ~ lagDv + successMaj + blockMaj + disappointMaj + success + block + disappoint + spendGapIMP + majControl + majInc + minControl + minInc + majorityQuality + minorityQuality + unopposed + majPresVote*presYear*majPresAgree + (1 | year) + (1 | district), data = data)
	summary(fourPart)

## this bit of code spins through all possible variable permutations ##	
## we also add an indicator for having been redistricted ##	

	core <- c('', ' + wins', ' + spendGapIMP', ' + majControl', ' + majInc', ' + minControl', ' + minInc', ' + majorityQuality', ' + minorityQuality',  ' + unopposed', ' + majPresVote', ' + presYear', ' + majPresAgree', ' + lagDv', ' + redist')
	append <- c('', ' + presYear*majPresAgree', ' + majPresVote*majPresAgree', ' + majPresVote*presYear',  ' + majPresVote*presYear*majPresAgree')
	
	length(core)
	length(append)	
	a <- c(0, 1)
	b <- c(0, 1)
	c <- c(0, 1)
	d <- c(0, 1)
	e <- c(0, 1)
	f <- c(0, 1)
	g <- c(0, 1)
	h <- c(0, 1)
	i <- c(0, 1)
	j <- c(0, 1)
	k <- c(0, 1)
	l <- c(0, 1)
	m <- c(0, 1)
	n <- c(0, 1)
	o <- c(0, 1)

	combs <- cbind(expand.grid(list(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o)))
	
	dim(combs) #should be 32768 x 15# 
	set <- list()
	for(i in 1:dim(combs)[1]){
		set[[i]] <- grep(1, combs[i, ])
	}

	tVals <- c()
	coef <- c()
	
	for(i in 1:32768){
		holdMain <- c()
		if(length(set[[i]]) > 0){
			for(c in 1:length(set[[i]])){
				holdMain <- paste(holdMain, core[set[[i]][c]])
			}
		}
		for(j in 1:5){
			mod <- paste('(majvote/100) ~ winsParty', holdMain, append[j], '+ (1 | year) + (1 | district)')
			model <- lmer(mod, data = data)
			coef <- c(coef, coef(summary(model))[2])
			tVals <- c(tVals, coef(summary(model))[2, 3])
		}
		cat('set', i, 'of 32768 run complete...\n')
	}

		
## now let's repeat with something we think works pretty well ##

	core <- c('', ' + wins', ' + winsParty', ' + majControl', ' + majInc', ' + minControl', ' + minInc', ' + majorityQuality', ' + minorityQuality',  ' + unopposed', ' + majPresVote', ' + presYear', ' + majPresAgree', ' + lagDv', ' + redist')
	append <- c('', ' + presYear*majPresAgree', ' + majPresVote*majPresAgree', ' + majPresVote*presYear',  ' + majPresVote*presYear*majPresAgree')
	
	length(core)
	length(append)
	
	a <- c(0, 1)
	b <- c(0, 1)
	c <- c(0, 1)
	d <- c(0, 1)
	e <- c(0, 1)
	f <- c(0, 1)
	g <- c(0, 1)
	h <- c(0, 1)
	i <- c(0, 1)
	j <- c(0, 1)
	k <- c(0, 1)
	l <- c(0, 1)
	m <- c(0, 1)
	n <- c(0, 1)
	o <- c(0, 1)

	combs <- cbind(expand.grid(list(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o)))
	
	dim(combs) #should be 32768 x 15# 
	set <- list()
	for(i in 1:dim(combs)[1]){
		set[[i]] <- grep(1, combs[i, ])
	}

	tVals1 <- c()
	coef1 <- c()
	
	for(i in 1:32768){
		holdMain <- c()
		if(length(set[[i]]) > 0){
			for(c in 1:length(set[[i]])){
				holdMain <- paste(holdMain, core[set[[i]][c]])
			}
		}
		for(j in 1:5){
			mod <- paste('(majvote/100) ~ spendGapIMP', holdMain, append[j], '+ (1 | year) + (1 | district)')
			model <- lmer(mod, data = data)
			coef1 <- c(coef1, coef(summary(model))[2])
			tVals1 <- c(tVals1, coef(summary(model))[2, 3])
		}
		cat('set', i, 'of 32768 1st run complete...\n')
	}

	hold <- data.frame(coef, coef1, tVals, tVals1)
	write.table(hold, file = 'allPermutationsResults.txt', sep = '|')

## now the picture ##
	
	pdf('figure1.pdf', width = 9, height = 9)
	par(mfrow = c(2, 2))

	plot(1, 1, type = 'n',
		xlim = c(range(0, coef1*100)),
		ylim = c(0, max(density(coef)$y)),
		xlab = 'Parameter Estimates',
		ylab = 'Density',
		main = 'Simplified Agenda Measure: Win Rate',
		axes = FALSE)
	axis(1)
	axis(2)
	polygon(density(coef), col = rgb(0, 0, 0, 0.25), border = rgb(0, 0, 0, 0.25))
	arrows(median(coef), max(density(coef)$y)/2, median(coef), 0, length = 0.125, lwd = 2, col = rgb(0, 0, 0, 0.75))
	text(median(coef), max(density(coef)$y)/1.75, paste('Median\n', round(median(coef), 3)), col = rgb(0, 0, 0, 0.75))

	plot(1, 1, type = 'n',
		xlim = c(range(0, coef1*100)),
		ylim = c(0, max(density(coef1*100)$y)),
		xlab = 'Parameter Estimates',
		ylab = 'Density',
		main = 'Spending Gap',
		axes = FALSE)
	axis(1)
	axis(2)
	polygon(density(coef1*100), col = rgb(0, 0, 0, 0.25), border = rgb(0, 0, 0, 0.25))
	arrows(median(coef1*100), max(density(coef1*100)$y)/2, median(coef1*100), 0, length = 0.125, lwd = 2, col = rgb(0, 0, 0, 0.75))
	text(median(coef1*100), max(density(coef1*100)$y)/1.75, paste('Median\n', round(median(coef1*100), 3)), col = rgb(0, 0, 0, 0.75))

	plot(1, 1, type = 'n',
		xlim = c(-5, 15),
		ylim = c(0, max(density(tVals)$y)),
		xlab = 'T-Values',
		ylab = 'Density',
		main = 'Simplified Agenda Measure: Win Rate',
		axes = FALSE)
	axis(1)
	axis(2)
	polygon(density(tVals), col = rgb(0, 0, 0, 0.25), border = rgb(0, 0, 0, 0.25))
	arrows(median(tVals), max(density(tVals)$y)/2, median(tVals), 0, length = 0.125, lwd = 2, col = rgb(0, 0, 0, 0.75))
	text(median(tVals), max(density(tVals)$y)/1.75, paste('Median\n', round(median(tVals), 3)), col = rgb(0, 0, 0, 0.75))
	polygon(c(1.65, 1.65, max(tVals1)+2, max(tVals1)+2), c(0, max(density(tVals)$y), max(density(tVals)$y), 0), col = rgb(0, 0, 0, 0.1), border = FALSE)
	arrows(1.65, max(density(tVals)$y), 15.5, max(density(tVals)$y), length = 0.125, lwd = 2, col = rgb(0, 0, 0, 0.5))
	arrows(1.65, 0, 15.5, 0, length = 0.125, lwd = 2, col = rgb(0, 0, 0, 0.5))
	text(5, 1.3, 
	paste('Prediction Region\nCoverage:', round(length(tVals[tVals > 1.66])/length(tVals), 3)), col = rgb(0, 0, 0, 0.75))

	plot(1, 1, type = 'n',
		xlim = c(-5, 15),
		ylim = c(0, max(density(tVals1)$y)),
		xlab = 'T-Values',
		ylab = 'Density',
		main = 'Spending Gap',
		axes = FALSE)
	axis(1)
	axis(2)
	polygon(density(tVals1), col = rgb(0, 0, 0, 0.25), border = rgb(0, 0, 0, 0.25))
	arrows(median(tVals1), max(density(tVals1)$y)/2, median(tVals1), 0, length = 0.125, lwd = 2, col = rgb(0, 0, 0, 0.75))
	text(median(tVals1), max(density(tVals1)$y)/1.75, paste('Median\n', round(median(tVals1), 3)), col = rgb(0, 0, 0, 0.75))
	polygon(c(1.65, 1.65, 17, 17), c(0, max(density(tVals1)$y), max(density(tVals1)$y), 0), col = rgb(0, 0, 0, 0.1), border = FALSE)
	arrows(1.65, max(density(tVals1)$y), 15.5, max(density(tVals1)$y), length = 0.125, lwd = 2, col = rgb(0, 0, 0, 0.5))
	arrows(1.65, 0, 15.5, 0, length = 0.125, lwd = 2, col = rgb(0, 0, 0, 0.5))
	text(5, 0.39, 
	paste('Prediction Region\nCoverage:', round(length(tVals1[tVals1 > 1.65])/length(tVals1), 3)), col = rgb(0, 0, 0, 0.75))
	dev.off()